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Abstract 

We introduce a new definition of distinguished trajectory that gener- 
alises the concepts of fixed point and periodic orbit to aperiodic dynamical 
systems. This new definition is valid for identifying distinguished trajecto- 
ries with hyperbolic and non-hyperbolic types of stability. The definition 
is implemented numerically and the procedure consist in determining a 
path of limit coordinates. It has been successfully applied to known ex- 
amples of distinguished trajectories. In the context of highly aperiodic 
realistic fiows our definition characterises distinguished trajectories in fi- 
nite time intervals, and states that outside these intervals trajectories are 
no longer distinguished. 
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This paper attempts to generalise the concepts of fixed point and 
periodic orbit to time dependent aperiodic dynamical systems. Fixed 
points and periodic orbits are keystones for describing solutions of 
autonomous and time periodic dynamical systems, as the stable and 
unstable manifolds of these hyperbolic objects form the basis of the 
geometrical template organising the description of the dynamical sys- 
tem. The mathematical theory of aperiodic dynamical systems is far 
from complete. In this context, this work deals with a general defini- 
tion that encompasses the concepts of fixed point and periodic orbit 
and which when applied to finite time and aperiodic dynamical sys- 
tems identifies special trajectories that play an organising role in the 
geometry of the fiow. 

1 Introduction 

In recent years the theory of dynamical systems has provided a useful framework 
for describing transport in fluid flows. Since the seminal work by Aref 1 on 
chaotic advection much progress has been made both in theory and applications. 
Dynamical systems techniques were first applied to Lagrangian transport in the 
context of two-dimensional, time-periodic flows [2] and stationary 3D flows such 
as the ABC flow 3 . More recently these techniques have been extended to 
describe aperiodic flows [H [7l |8] and flnite time-dependent flows, such as those 
rising in geophysical applications I12j . However, the mathematical theory 
for both aperiodic time-dependent flows and finite time aperiodic flows is far 
from being completely developed. 

For stationary flows the idea of fixed point is a key for describing geomet- 
rically the solutions. Fixed points may be classified as hyperbolic or non- 
hyperbolic depending on their stability properties. Stable and unstable man- 
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ifolds of hyperbolic fixed points organise the phase portraits of the flow away 
from the region close to the fixed points [5J |2] ■ These manifolds comprise re- 
spectively the trajectories that approach the fixed points as time tends to plus 
or minus infinite. As they are formed of trajectories they act as barriers to 
transport as particles cannot cross them without violating the uniqueness of 
the solution. They are useful because they allow qualitative predictions for 
the evolution of sets of initial conditions avoiding explicit integration of initial 
conditions on the whole domain. Hyperbolic fixed points and their stable and 
unstable manifolds are the basic notions used for the geometrical description of 
flows in autonomous dynamical system. 

The concept of fixed point is extended to time periodic flows by means of 
the Poincare map, as periodic orbits with period T become fixed points of the 
Poincare map. For hyperbolic periodic orbits there exist also stable and unstable 
manifolds that are geometric objects that organise the global dynamics. Again 
they are respectively the sets of orbits asymptotically approaching the periodic 
orbit as time tends to plus or minus infinity. 

Aperiodic fiows are still poorly understood, as theory that is well established 
for autonomous or periodic fiows does not apply to them directly. For instance 
there exists efforts in the mathematical community to extend the well known 
concept of bifurcation for stationary flows to non- autonomous systems [HI llOj . 
To gain insight on the geometrical structure of aperiodic flows, concepts such 
as Lyapunov exponents are used, however these are defined strictly on infinite 
time systems. Realistic fiows, like those arising in geophysics or oceanography, 
are not infinite time systems and for their description, finite time versions of 
the definition of Lyapunov exponents such as Finite Size Lyapunov Exponents 
(FSLE) [14 and Finite Time Lyapunov Exponents (FTLE) jUl [16] are used. 
Special trajectories, such as detachment and reattachment points jl9j . are ob- 
served in highly aperiodic or turbulent fiows. In particular these separation 
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trajectories occur on the boundaries in simplified ocean models |21j and also in 
technological applications in air foil design ^U\. Recent articles by Ide et al and 
Ju et al [T71 [T5] referring to these special trajectories introduce the concept of 
Distinguished Hyperbolic Trajectory (DHT) which encompases not only trajec- 
tories on the boundaries but also special trajectories in the interior of the flow. 
DHT are hyperbolic trajectories that, like hyperbolic fixed points and periodic 
orbits, have stable and unstable manifolds that are key for describing geomet- 
rically the solutions on the phase space. This generalisation is an important 
step- forwards in the study of aperiodic fiows, as it is a powerful tool for describ- 
ing transport in realistic oceanographic fiows [TTJ |T5J |T31 HZj. Distinguished 
hyperbolic trajectories as defined in [T71 [T5] are computed from hyperbolic in- 
stantaneous stagnation points (ISPs) by means of an iterative procedure. If 
instantaneous stagnation points bifurcate and do not persist for all times the 
technique developed in |17lll8j cannot be applied in those time intervals, leaving 
many questions unanswered, such as what happens to the distinguished trajec- 
tories at those times, for distinguished hyperbolic trajectories are trajectories, 
and as trajectories exist at all times. In fact Ref. pLl, provides examples of 
vector fields with exact distinguished hyperbolic trajectories that exist on time 
intervals without hyperbolic ISP. Refs. [TTJ [T^] discuss the impossibility of this 
technique for tracking DHTs after ISP bifurcations and as a consequence the 
difficulties in establishing whether DHTs obtained at different times are part of 
the same trajectory or not. 

In this paper, following ideas discussed in [Til El Hill we propose a new 
definition of Distinguished Trajectory (DT) which generalises the concepts of 
fixed point and periodic orbit to aperiodic fiows. We have taken the liberty of 
calling them Distinguished as in |11|, I17|, I18j , since although the definitions are 
not strictly equivalent, it is found that the studied hyperbolic trajectories are 
encompassed by both definitions. We remark that our notion has the advantage 
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over the method proposed in [T71 [TH] that the DTs may be computed without 
the presence of hyperbolic instantaneous stagnation points. Our definition does 
not depend on the dimension n of the space on which the vector field is defined 
and is valid both for hyperbolic and non-hyperbolic types of stabilities. Non- 
hyperbolic DTs have not been studied in [T71[TH], and in this sense our definition 
is broader than that proposed there. In particular, we will show that exact non- 
hyperbolic periodic orbits fall within the category of distinguished trajectories. 
Trajectories of this type could be of special interest for their applications in 
oceanography, as they are related to eddies and vortices. Ocean eddies are 
well studied |28j . Frequently they are long lived, and water trapped inside can 
maintain its biogeochemical properties for long time, being transported with the 
vortex. In steady horizontal velocity fields, the presence of closed streamlines is 
the mathematical reason for the isolation of the vortex core from the exterior 
fluid. In two-dimensional, incompressible, time-periodic velocity fields the KAM 
tori enclose the core, a region of bounded fluid particle motions that do not mix 
with the surrounding region 4 . But how to deflne an eddy from the Lagrangian 
point of view in aperiodic flows? This is still an open question |27l 129] for which 
we will discuss new possibilities suggested by the deflnitions given in this article. 

The structure of the paper is as follows. Section 2 introduces the definition 
of distinguished trajectory and explains its motivation in the context of ID 
examples. Section 3 explains the algorithm used to verify the applicability of our 
definition of distinguished trajectories to the solutions of the periodically forced 
Duffing equation. Details about technical issues arising from implementation 
of the definition are given. Section 4 reports the results obtained in several 
other 2D and 3D examples, both periodic and non-periodic, hyperbolic and 
non- hyperbolic. Section 5 discusses results on realistic fiows. Attention is paid 
to open questions on distinguished trajectories such as those mentioned above 
and pointed out in Refs. [121 |TT]. Finally, section 6 presents the conclusions. 
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2 Distinguished trajectories: a definition 

We start by recalling the definition of distinguished hyperbolic trajectory pro- 
vided in [17'. Given the system: 



Let x{t) be a trajectory of Eq. that remains in a bounded region for all 
time. Then x(t) is said to be a distinguished hyperbolic trajectory if: 

1. it is hyperbolic, 

2. there exists a neighbourhood B in the flow domain having the property 
that the DHT remains in B for all time, and all other trajectories starting in B 
leave B in finite time, as time evolves in either a positive or negative sense, 

3. it is not a hyperbolic trajectory contained in the chaotic invariant set 
created by the intersection of the stable and unstable manifolds of another 
hyperbolic trajectory. 

Remark 1 // the data spans only a finite time interval, then the DHT cannot 
be determined uniquely. Instead, there is a small region in B where the DHT 
can exist. 

In [17] this setup is extended to general vector fields as follows. Coordinate 
transformations are sought which put the system in the form of Eq. 1 and then 
the previous definition is applied. 

We give now our definition of distinguished trajectory for a general vector 



dx 



Dx + 5^'^(x,t) xe W 



,71 



(1) 



field: 



v(x,t), X e M", t e M 



(2) 
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We assume that v(x, t) is C" (r > 1) in x and continuous in t. This will allow for 
unique solutions to exist, and also permit linearization, although linearization 
will not be used in our construction. 

Before giving our definition of DT, we first need to introduce some notation 
and to make some definitions. Let x{t) denote a trajectory of the system Q 
and denote its components in M" by (xi, 2:2, For any initial condition x* 

in an open set B C M" , consider the function M : S ^ M 




M is the function that associates to each initial condition x* in B the arclength of 
the trajectory that passes through x* at time t* . The arclength of the trajectory 
is considered over its projection in the phase space (xi, X2, •■•x„) and depends on 
t* and r. As the function M is defined over an open set it does not necessarily 
attain a minimum, but if it does, the minimum is denoted by min(Af (x*)t* t-). 

Definition 2 (t- Distinguished trajectory). A trajectory j{t) of Eq. ^ is r- 
distinguished at time t* if there exists an open set B around j{t*) on which the 
defined function M{x*)f,T has a minimum and 

min(Af(x*)t.,,) = M(7(r))t.,,. (4) 



2.1 A discussion of the definition 

The elements of the above definition deserve a detailed justification. We il- 
lustrate our explanations with examples in ID. First we consider an example 
taken from |17l 122] . It is the linear one-dimensional non-autonomous dynamical 
system given by: 

dx 

- = -x + t. (5) 



For this example we consider the DHT reported in |T7], which is given by 2: = 
i — 1 . This is the particular solution of the linear equation ([s]) towards which 
all trajectories decay. The solution through the point x* at i = is given by, 

x{t) ^t-l + e-'{x* + 1). (6) 

Figure [TJl) displays several trajectories starting at times ranging from t = to 
t — A and figure ^p) displays the same but starting at time t = —4 (note that 
in this case part of the trajectories are out of the displayed domain). For each 
initial condition the function M provides the length of the projection of the 
trajectory over the x-axis in the range of times [— r, r]. Geometrically it is clear 
that in this example the function M should have a minimum for a certain x value 
and that this value depends on r. Ideally the minimum of M should coincide 
with the position of the DHT at t = 0, however this would not be possible if 
in the definition of M only positive times were considered, i.e. if the limits 
of the integration were (0,r) the dashed trajectory in figure [iji) would have a 
lower projection in positive times than the particular solution. An analogous 
problem would be encountered were only negative times considered, that is if 
the limits of the integration would have been (— t, 0). To determine precisely 
the position of the DHT at t — 0, both positive and negative times must be 
considered in the definition of M. Figure lb) confirms that with this choice 
the dashed trajectory cannot be distinguished, as it increases its projection in 
negative times. Figure [2^) displays the function M(a;*)t=o,r evaluated along 
the trajectories for several t values. Figure |2]d) displays the position of the 
minimum of the function Mt=o,r as a function of r. These minima correspond 
to the positions of the r-distinguished trajectories at t = and as r increases 
they approximate the coordinate of the DHT at this time, which is at x* = — 1. 
The pair (t',x') formed by the time at which AI is computed and the value of 
the coordinate x' to which the minimum of the function Mfi j. converges for 
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increasing t is called the limit coordinates. The graphic [2p) illustrates the idea 
of approaching a point (io,xo) of the distinguished trajectory by means of the 
limit coordinates. In practice the convergence to the limit coordinates cannot be 
examined in the limit r ^ oo, either because it is impracticable in a numerical 
implementation, or because in the large r limit errors accumulate, or simply 
because the dynamical system is defined by a finite time data set. For these 
reasons the convergence to the limit coordinates will be tested up to a finite r. 

Figure ^p) raises the question: what controls the rate of the convergence 
of the minima of M to the coordinates of the DHT? It is hard to answer this 
question rigorously for a vector field as general as in Eq. (|2|. However, some 
insight may be provided by particular examples. For instance the system 



has the same DHT as ([5]). Its solution through the point x* at i = is given by 



Here the decay of the solution towards the DHT is faster due to the presence 
of the exponential term e"^*. Figure [s] shows that in this case the rate of the 
convergence of the minima of M towards the coordinates of the DHT at time 
t = is also faster than before. However there exist systems in which the 
exponential decay of the solution is not a determining factor affecting the rate 
of the convergence of the minima of AI to the coordinates of the DHT. For 
instance, in autonomous systems fixed points are the DTs, and clearly they are 
minimizers of M for any t > whatever is the exponential rate of growth or 
decay of the nearby solution. 

In these examples the function M has a unique minimum, but as we will 
see the situation will not always be so simple when nonlinearities are involved 
in the vector field. Also it is important to notice that the function M obtained 



dx 



= -2x + 2t-l, 



(7) 



x{t) =t-l + e-^*{x* + 1). 



(8) 
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at different r values has been used to obtain the limit coordinates (^qjXq) and 
that these approach the xq coordinate of DHT at a given time to (here to = tg). 
Once this is obtained, approaching the DHT at later times = to + kAt would 
require applying the same procedure to get the limit coordinates (tj^jxjj,). We 
remark here that the proposed algorithm does not ensure that the set of limit 
coordinates (i'j,,x'j,) are in fact part of a trajectory. Later we will see that in 
practice, in many examples these points approach a true trajectory, however in 
realistic aperiodic flows this has to be verified a posteriori. These considerations 
lead us to the definition of a Distinguished Trajectory. 

Definition 3 (Distinguished trajectory). A trajectory ^(t) is said to be Dis- 
tinguished with accuracy e (0 < e ) in a time interval [tQ,tj\j] if there exists a 
continuous path of limit coordinates (<',x') where t^ E [^Oi^Jv]; such that, 



In the numerical exploration of this definition we will replace the continuous 
path of limit coordinates (t',x') and the continuous trajectory j{t) by discrete 
representations (t[.,x'j,) and 7(t[.) where to < < tN- By definition [s] any 
trajectory is distinguished for sufficiently large e, however the interesting distin- 
guished trajectories are those for which e is close to zero, which means it is of 
the order of the accuracy in which 7(<5j) and x'(t5,) are numerically determined, 
or zero, if an exact expression is known for both. 

Underlying definitions 2 and 3 is the geometrical idea that distinguished 
trajectories, which act as organising centres of the fiow in phase space, are 
those that "move less" (in a certain sense) than other nearby trajectories. This 



(9) 



Here \ \ ■ \ \ represents the distance defined by 
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property of "moving less" is satisfied by minima of the function M as it measures 
the length of the displacement in phase space of a trajectory forwards and 
backwards in time. In fact this property is related somehow to property (2) of 
the definition provided in [17 and presented at the beginning of Section 2, as 
the trajectory that "moves least" is not expected to leave the neighbourhood B. 

Definitions[2]and[3]are made for a general dynamical system in any dimension 
n. The purpose of this paper is the exploration of these definitions, but more in 
an illustrative than demonstrative way, as it is impossible to provide examples 
for every possible n, and one cannot deal with every possible example at a given 
n. Even if one wants to provide a rigorous formal proof that the definition 
recovers specific trajectories such as periodic orbits (it is not obvious that in 
general they have to satisfy our definition) , this has to be done with some further 
hypotheses on the vector field and proofs will not be valid beyond the assumed 
hypotheses. Therefore we restrict the discussion to dimensions up to 3, as these 
are the dimensions important for geophysical fiows, which are what originally 
motivated the definition. However it is sensible to make the same definition for 
any dimension n, as it is clear that it works for autonomous systems of any 
dimension. Fixed points are the kind of trajectory expected to be recovered by 
the definition and they do not move at all in the phase space. For these Af — 0, 
while M > for any other trajectory in the neighbourhood which is not a fixed 
point. 

We conclude this section with some remarks. First, it is not guaranteed 
a priori that for an arbitrary vector field, satisfying only some rather general 
hypotheses such as those of Eq. (|2]), the function M will have a minimum, how- 
ever this is not a problem from the point of view of the definition. For instance 
the same thing happens for general non-linear autonomous systems. In these 
systems fixed points are perfectly defined although one does not know a priori if 
such points exist for arbitrary examples. If they exist, it is possible to find them 
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by either solving the nonhnear equation v(x) = 0, or by applying definitions |2] 
and[3j In the same way one does not know a priori if distinguished trajectories 
exist for a general vector field however if they exist they can been found with 
the tools proposed in this article. Second, even if a path of limit coordinates is 
found it is not guaranteed that it will be a trajectory, although if that is the 
case then from definition [3] follows that this trajectory is distinguished. Third, 
one might think that if limit coordinates are found at to that approach with 
great accuracy a point of an existing DT, then the iterative procedure described 
above for finding a set of limit coordinates (<5,,,x5,) approaching the DT at later 
times is an unnecessary computational effort, as those coordinates could have 
been equally well obtained by integrating forwards the initial data. However 
there exist examples such as a hyperbolic DT in dimension greater than one 
with at least ID unstable manifold, that cannot be integrated like this, as the 
integrated trajectory will eventually leave the neighbourhood of the DT through 
the unstable manifold no matter how small the initial error is. In summary the 
proposed methodology based on limit coordinates provides a systematic way of 
finding DT, which can be elusive and difficult to obtain. We will discuss these 
issues in detail in later sections. 

3 A numerical algorithm 

In this section we propose an algorithm for computing a path of limit coordinates 
in a time interval, and we verify that it is close to a DT of a known example. For 
this purpose we calculate, at increasing r values, the minimum of the function 
Mt=o,r(x) for X in an open set in M". The method is illustrated in a 2D case, 
the periodically forced Duffing equation 

i = y, 

y = X ~ + e sm{t) , (10) 
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where e is a small parameter. The hyperbolic fixed point of the unperturbed 
autonomous system (i.e., e = 0) is at the origin x = (0,0). For small e, it is 
possible to compute by perturbation theory (see [23]), the following periodic 
trajectory which stays close to the origin: 



e/smt\ £3 /2sm^<+ fsmtcos^ A 

XD//T(i) = -- , o 9 +£>(&)■ (11) 

"^""'^^ 2\costJ 40 Vicosn + Ssin^tcosiy ^ ' ^ ' 



For e = 0.1, Eq. (11) is accurate up to the fifth digit. This trajectory is 
identified as Distinguished in Ref. 17J, for this reason we have labelled it a 
DHT. Substituting the expression, 

x= (a;,y) =Xi5HT(t) + (ei,6) (12) 



into Eq. (10 1 and by dropping the nonlinear terms one finds that the linearized 
equations have two linearly independent solutions in terms of which the time 
evolution of the components (^1,^2) is: 



(ei,6) = «e*( \^/"^)+^^-'[^\^/"^ ]+0{e'). (13) 



Eq. (13) confirms the hyperbolicity of the solution ( [Tl| . 



This explicit expression for the distinguished hyperbolic trajectory is a bench- 
mark for testing the utility of our definition. The procedure starts by determin- 
ing the coordinates of yiuHT at time t = 0. We consider the open set I? C M^, 
defined by P = (-0.2, 0.2) x (-0.2, 0.2) and in the function 7\fi^o,r(x) we take r 
to be 2. Figure |4] displays a contour plot of Aft^o,T=2(x) which has a minimum 
at X = (0, —5.7057 • 10^^). Aft=o,T=2(x) quantifies displacements of particles in 
phase space, and its minimum corresponds to the initial condition that "moves 
less" over the r interval [—2, 2]. As noted in the previous section, when the value 
of r is increased, the position of the minimum gets closer and closer to the coor- 
dinates of the DHT at t = 0. Figure [s] shows contour plots of the function M for 
several t values. Fig. [sji) displays a typical hyperbolic structure for M for r = 5 
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where the directions of the stable and unstable manifolds are easily recognised. 
In Fig. [5^) the function M has a unique minimum at x = (0,-4.979 • 10"^) 
while in Fig. [sja) there appear several minima for r = 10. The global minimum 
in this picture corresponds to x = (0,-5.0042565261 • 10^^). Figure [6^) com- 
pares the a:-coordinate of TioHT as a function of time with trajectories having 
initial conditions at the global minima of il/f=o,T=2 and Mt=o,r=io- Taking as 
initial condition the global minimum of Aft=o,r for r = 10 provides a trajectory 
that stays close to TioHT for a longer time interval than for r = 2, which con- 
firms that larger t- values more closely approach the coordinates of the DHT. 
Fig. [5j:) displays the contour plot of A/(=o.t-=5o(x). Its global minimum is at 
X = (0, -5.0037606418 • 10"^)^ xhe associated trajectory depicted in Fig. ^) 
shows that this initial condition tracks the DHT for a longer time interval than 
those obtained for t = 2 and 10, however the figure shows that the integration 
of the DHT in (—50,50) is not possible. In fact the associated trajectory stays 
close to the DHT only in the time interval (—20, 20). This confirms that results 
obtained for t = 50 are the same as those obtained for r = 20. In practice for a 
finite precision numerical scheme, such as a 5th order Runge Kutta used here, 
the approach to the DHT has an upper bound depending on r. This occurs 
because the stable and unstable manifolds of the hyperbolic trajectory magnify 
any initial error in either negative or positive time and beyond this r-limit nu- 
merical errors dominate. The convergence towards the DHT is confirmed in Fig. 
[7] which displays the evolution of the coordinates x and y of the global minimum 
of M as a function of the parameter r. 

New minima appearing in Figs, [sj)) and c) relate to the existence of different 
T-distinguished trajectories. As illustrated in Fig. |6|d), they correspond to 
trajectories which stay close to TioHT in a small time range contained in the 
interval — r < t < t, but which later fly apart from the DHT. 

We now describe a numerical scheme to compute a path of limit coordinates. 
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The algorithm has the following steps: 

1. Step 1. Discretize the domain T) at the initial time t = a,t which one 
wishes to compute a DT. For instance, the grid size of this domain in Fig. 
|4]is 101 X 101. The function M is evaluated at each grid point for a given 
To- 

2. Step 2. Search for the local minima of Mto.n, in the interior of the grid. 
These minima approach the coordinates of To-distinguished trajectories 
within the accuracy of the grid. In what follows we restrict our description 
to the case of a unique minimum, as this simplifies the description; the 
procedure is easily generalised to the case of multiple minima. 

3. Step 3. Improve the approach of the coordinates of the ro-distinguished 
trajectory up to precision S. For this purpose build up a 3" grid centred 
on the candidate point provided by step 2, (for the 2D case this is a 3 x 3 
grid as Fig. [s] illustrates) , setting the distance between nodes equal to S. 
Then evaluate Mtg^ro at the points of the (5-grid. If the minimum of Mtg,To 
is in the interior of the grid, then the coordinates of the rg-distinguished 
trajectory are known to within 6 accuracy. Otherwise the (5-grid must 
be rebuilt centred on the boundary point where the minimum has been 
located, and Mj„ ^-o must be re-evaluated in the new (5-grid. This procedure 
stops when the minimum of Mt^ T^ is in the interior of the grid. 

4. Step 4. Computing the limit coordinates at time to- Define a sequence of 
increasing r-values as follows: ti = tq + Ar and T2 = tq + 2Ar. Then 
evaluate Mtg^ro^ -^^to,Ti and Mtg,T2 on the (5-grid. If the minimum is at an 
interior position for the three cases, then we consider that limit coordi- 
nates have been found within S accuracy. We note that this is a necessary 
but not sufficient condition as one does not know a priori the conver- 
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gence rate to the distinguished trajectory. Although this criterion could 
be strengthened, it has been tested and found to be adequated for the 
examples explained in subsequent sections. If the condition defined above 
of having a minimum at an interior position for the sequence of t- values 
is not satisfied then, after replacing tq by ri, we return to step 3 and then 
to step 4. The loop between steps 3 and 4 is stopped when the condition 
of step 4 is satisfied for some . 

5. Step 5. Compute the limit coordinates at time ti = tg + At. Once the limit 
coordinates have been approached at time to- they are integrated forward 
numerically up to time ti. If the limit coordinates converge to a hyperbolic 
DT with an unstable manifold, the position x(ti) obtained should deviate 
from the position of the DT at time ti. In order to correct this, the pro- 
cedure described above is repeated from step 3 onwards. For that purpose 
in the definition of M, to is replaced by ti and the r-value is reset to tq. 
The coordinates x(ti) are the first approximation to the ro-distinguished 
trajectory at time ti. Once the limit coordinates are found for time ti 
it is possible to repeat the procedure to locate them at successive times 
t2 , , . . . , tjv . 

The algorithm requires as inputs: an explicit expression for the dynamical 
system Q; the definition of the domain V C M"; the initial and final times to, 
t^ at which DTs are required, and the time step At for intermediate times; the 
initial tq and the increment At; the precision d. As an output the algorithm 
gives a path of limit coordinates at the selected times tk- 

Next we discuss in more detail some technical aspects related to the imple- 
mentation of the above algorithm. Steps 1 and 3 require evaluating Mt„,To ^ 
defined in Eq. (|3|. We explain how this is done for the contour plots displayed 
in Figs. [4] and |5] which refer to the system (10) at to = 0. Fig. [9] shows 
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a schematic projection onto the plane of a possible trajectory x(t) of the 
system from — r to r. As it was obtained numerically, only a finite number of 
points (L) appear. This picture suggests the following discrete version of Eq. 
([3| for M: 

where the functions xj (p) and yj (p) represent a curve interpolation parametrized 
by p, and the integral 



/ dxj{p)Y , ( dyj{p)^ ^ 



Pi 



, , . , , I dp (15) 

dp J \ dp J 

is computed numerically. In our case we use the Romberges method (see f25j) 
of order 2K with K = 5. It is clear that the accuracy of the evaluation of M 
will depend on the number of points on the trajectory L, which is controlled by 
the size of the time step, h, of the integrator (a 5th order Runge Kutta method) 
and on the interpolation scheme between points. Two interpolation methods 
are compared in tables ([l]) and Results in table (jlj are obtained with linear 
interpolation between nodes. Results in table ^ correspond to the interpolation 
method used by Dritschel 26J in the context of contour dynamics, which has 
been successfully applied in |23j to the computation of invariant manifolds for 
aperiodic flows. This method interpolates a piece of the curve in Fig. |9]between 
consecutive nodes as follows: 

(p) = Xj + ptj + rjj {p)nj (16) 

for Pi = 0<p<Pf = l with Xj(0) = Xj and Xj(l) = Xj+i, where: 

tj = (aj,&j)-xj+i-x„ t,eR^ (17) 
Uj = {-bj.aj), rij e (18) 
Vj{p) -= PjP + f3jP^ + ljp\ (19) 
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The cubic interpolation coefficients fij, (3j and 7^ are 

fij — —-djKj — -djKj+i, (3j = '^djKj, 7j = -dj{Kj+i — Kj), 
where dj = — Xj| and 

is the local curvature defined by the circle through the three points, 
and Xj^i. 

Tables ([l]) and ^ show the errors in the computed lengths of the ellipses 
for different ratios of major to minor axis. The reference length is that obtained 
with GNU Octave version 2.1.73, as it provides 16 correct digits for the known 
circumference. The tables confirm that the Dritschel's method is superior to 
linear interpolation and it is the one used to compute the function M. In the 
trajectory from — r to r the number of points L is determined by the time step 
size of the Runge Kutta method which is set to 10"^. 

Another important element of the algorithm needing discussion is the value 
of the input parameters, in particular of tq and At. It is clear from Fig. [Sjthat 
large r values are not convenient as they increase the roughness of the function 
M and several local minima may appear in the neighbourhood of a DHT that 
correspond to trajectories that stay close to it for some time. On the other hand 
it is clear that sufficiently large r values are required to fix the coordinates of 
the DHT to within prescribed accuracy. Combining these observations suggests 
the use of relatively small values for the initial tq. In the example above tq = 2, 
provides, as a starting point, a smooth M as that of Fig. |4] The increments 
should not be large. In practice we have chosen At = to/2. This prevents 
from stepping to a too rough Af before getting close enough to the sought after 
DHT. Some of the local minima appearing in Fig. ^p) are just apparent and 
disappear with a more refined grid. However, as already observed, others belong 
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to true T-distinguished trajectories, which are secondary and can be avoided if 
the increment of the r-values is conveniently small. These choices are found to 



be appropriate for determining with great accuracy the DHT in ( 1 1 1 by means of 



a path of limit coordinates. Figure 10 1) represents both the analytical DHT and 



the numerical limit coordinates and Figure 10 d) displays the distance between 
the exact and the numerical approach, confirming that the DHT in (11 1 is also 
a DT in the sense of our definition [s] with accuracy e = 10^^. Other parameters 
in the algorithm are: S — 10^^, step size in the Runge Kutta method h = 10~^, 
io = 0, = 6, and At = 0.01. To locate the DHT with accuracy 5 = Ur^ 
requires increasing values of r up to 15, which is near the limit of the integration 
method. Figure [TT] shows the maximum r required at each t^. 

4 Applications to exact examples 

In this section we apply the algorithm explained in the previous section to 
selected examples. 

4.1 A non- hyperbolic distinguished trajectory 



The unperturbed autonomous system ( 10 1 obtained with e = has non- hyperbolic 



fixed points at (—1,0) and (1,0). Obviously these fixed points correspond to 
DTs which are also r-distinguished trajectories for all r > 0. For the period- 



ically forced system (10) with small e it is possible using perturbation theory 
to find periodic solutions close to these fixed points in a manner similar to the 
analysis of the hyperbolic example made in the previous section. For instance 
close to the point (1,0) we find the periodic trajectory: 

...T(t) = -('\ +e(''^'\ +3s^( ^.7'^ Uo(e3). (20) 
\0/ \costJ \— smtcostj 

This solution has not been considered Distinguished in previous works |17Lll8j . as 
these have been focused on hyperbolic trajectories and this solution, as is proved 
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next, is not hyperbolic. However, in anticipation of its having the Distinguished 
property, we have labelled it DET for two reasons. One is that it is periodic, 
and we expect periodic orbits to be distinguished, and second is that it is in 
clear correspondence to the elliptic fixed point (-1,0) in the case e = 0, and fixed 
points are DTs. 



To determine the stability of (20) we proceed as before, by substituting the 
expression 

X = (x, y) = ^DET{t) + (6,6) (21) 



into Eq. (10 1. We find that the linearized system at order is 

dt 
d^2 



dt - ^''^ 



dt - ^''^ 
Therefore the linearized fiow around TiDErit) evolves according to: 

which clearly is not hyperbolic. Here a and a* are complex conjugate numbers. 



We apply our algorithm to determine the limit coordinates approaching (20 1, 
as we want to verify whether definition [3] also works for time-dependent non- 
hyperbolic solutions. The following input is considered: V — (—1.2,-0.8) x 
(-0.2,0.2), To = 2, At ^ 1, S = 10"*, to = 0, = 6, and time step 10"^ foj. 
the Runge-Kutta integrator. We note that the accuracy 6 is not as demanding 
as before, since now the exact TioET for e — 0.1 is only accurate up to the 
third digit. Figure [12] shows a rather different structure for the function M. An 
important feature is the smoothness of M close to the DET even for large r. 



In figure 12 d) the differences between the rather flat region around the position 



of the DT given by (20 1, which appears in the dark tone, and the roughness 
of the outer part are remarkable. The irregularity of this region suggests that 
inside it nearby trajectories follow rather different paths as happens for chaotic 
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motions, while the regularity of the central core suggests the existence of trapped 
trajectories circling around the DET. From this perspective the function M for 
large r seems a useful tool for fixing the boundaries of a Lagrangian eddy, 
different to the methods proposed in [271 121] • 

Figure [13] shows the rate of convergence to the global minimum of AI in 
the domain I? as a function of r. The convergence towards the coordinates of 
the DT is oscillatory and rather slow since t values up to 600 are required. A 
slight difference between the exact coordinates of the DT and the numerically 
computed limit coordinates is evident, however we note that these differences 
are consistent with the precision to which the exact DT is known, which is only 



to the third digit. Figure 14 and more specifically figure 14 d), confirms that 



the exact expression in Eq. ( 20 1 is in fact a distinguished trajectory according 

1-3 



to our definition [3] with accuracy e = 4 • 10 



Figure 15 1) shows a forward and backward integration along the time interval 
(^50,50) taking as initial data the limit coordinates supplied by our algorithm 
at time — Q, and compares it with the exact solution of the DT. From figure 



15 d) it can be seen that this trajectory evolves close to the exact solution in 



the entire time range. This result shows that contrary to what happens near 
hyperbolic trajectories, near non hyperbolic trajectories, small error does not 
amplify and as consequence, once a DT is known to exist it could have been 
computed simply by integrating forwards and backwards the limit coordinates 
found at a given time ij.. However one needs to be careful here, as a trajectory 
is not necessarily distinguished at all times, and for it to be properly called 
distinguished, it should be verified that it stays close to the limit coordinates 
in the whole time interval, and therefore one cannot avoid computing limit 
coordinates along the time interval in this case either. We will return to this 
point in the next section. 
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4.2 The rotating DufRng equation 

Next we analyse the aperiodic hyperbolic distinguished trajectory of a system 
already studied in [23 , the rotating DufRng equation: 

fji \ _ f sin 2tJt cos 2ujt + uj \ /' 771 
772 I \ cos 2ujt — oj — sin 2ujt I \ rj2 

+ {e sint- [cos ujtTji- sin Ljti]2f )( l^^J^l ). (25) 



cos ujt 

This DufRng equation is quasi-periodic in time when the rotation rate lu is 



irrational. It is obtained from the system (10 1 by applying the rotation x = 
R{t)T], where 

i?(t)=('^°^'^! -sina;i\ 
^ ' y sm Lut cos ivt J ^ ' 

The DHT can also be obtained through the coordinate transformation: 

- mr^^DHTit). (27) 



Figure 16 and in particular Fig. [I6p), confirms that the DHT (27 1 is also a 
DHT according to our definition [sjwith accuracy e = 4 • 10^^. 

4.3 A 3D extension of the DufRng equation 

In this section we apply our definitions to an example in higher dimension. In 
particular we consider a 3D extension of the Duffing equation: 

y = a; - + esin(i), (28) 
i — z + £sin(t). 
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The hyperbolic fixed point of the unperturbed autonomous system (i.e., £ = 0) 
is at the origin x = (0, 0, 0). The solution for small e becomes: 

^DHT{t) = 




e 
40 



3/2 sin^ i + I sin t cos^ t 



fcos^i + Ssin^tcost I +0{e^). (29) 




The numerical scheme explained in Section 3 is easily adapted to higher 
dimensions. However some changes must be made. The computation of M re- 
quires approximating lengths of trajectories which in 3D needs an interpolation 



scheme different to that of Eq. (16 1, which is only valid in Mr. We consider 
the linear interpolation instead. This interpolation evaluates the function M 
satisfactorily if trajectories are represented by a large number of points. This is 



achieved by using a Runge-Kutta method with time step h = 10^^. Figure 17 
indicates the evolution of coordinates associated with the minimum of M as a 
function of t (solid line). The dashed line corresponds to the exact perturba- 
tive solution. There is evident a clear convergence towards the exact position 
although there is a significant jump in the asymptotic behaviour beyond r ^ 50. 
This jump is due to round off errors in the determination of M for large r. The 



third equation in (28 1 is just a linear equation and for this reason solutions which 
are in the neighbourhood of the DHT have z-coordinate growing exponentially 
in backwards time. Thus for large t values, the evaluation of M is made along 
very long trajectories in the z-coordinate, which are underrepresented by points 
sampled every h = 10^^ (see table I) and where lengths are badly calculated 
by adding up very small and very large (and inaccurate) numbers. In spite of 
this, figure [18] confirms that the exact distinguished trajectory can be accurately 
obtained with our methodology and that for t < 50 errors are within the ex- 
pected margin. The remaining input parameters used in Figs. [17] and [18] are: 
V = (-0.2,0.2) X (-0.2,0.2) x (-0.2,0.2), To = 2, At = 1, (5 = IQ-^, step size 
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h = 10~^ in the Rungc Kutta method, to = 0, = 6, and At = 10~^. 

As we explain next, the computational demands made by this example are 
considerably larger than they were for the previously considered 2D examples. 
When determining a DT, most of the CPU time is spent computing the value of 
M on the (5-grid displayed in Fig. [8j The number of neighbours of the interior 
point grows with the dimension n as 3", therefore when the problem increases 
its dimension from n to n + I, the computational demands are multiplied by 3. 
Another factor that contributes to increased computational time is the decrease 
of the Runge Kutta time step h in the evaluation of trajectories on the (5-grid. 
This increases the number of points in the trajectory (and therefore the number 
of operations) with respect to the previous Dristchel approach by a factor 100. 
This factor is partially balanced by the fact that for the same number of points 
the arclength is computed more rapidly with the linear than with the Dristchel 
interpolation. 

5 Application to vector fields defined as finite 
time data sets 

In this section we explore definition [3] for a highly aperiodic 2D flow in which 
the vector field is defined as a finite time data set. In particular we consider the 
output of a quasigeostrophic wind-driven double gyre model in a regime already 



studied in [TTJ [T^] . Details of this model may be found in jT^l [H] . Fig. 19 shows 
a typical output for the streamfunction provided by this model. The velocity 
data set is obtained on a 1000 km x 2000 km rectangular domain and spans 
4000 days. This interval is considered for a fluid started from rest and allowed 
to spin for 25000 days. Free slip conditions are considered for the velocities 
on the boundaries and the wind stress curl is 0.32 dyn/cm^. The equations of 
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motion for this system are given by: 

X = Vx{x,y,t) ^ (30) 
dy 

y = ""vi^^V^^) = (31) 

and the variables a; and y are in the rescaled domain [0,1] x [0,2]. Here the 
velocity fields and Vy are provided as a finite time data set and are interpo- 
lated using bicubic interpolation in space and 3rd order Lagrange polynomials 
in time. This method has been reported to be good enough for integrating tra- 
jectories in [2|]. We will focus our analysis in the time interval [0,900] in the 



area marked by a rectangle in Fig. 19 for which jl2j reports the computations 



of several DHTs. In [T^] distinguished trajectories are computed by means of an 
iterative algorithm which is initialized on a hyperbolic instantaneous stagnation 
point (ISP). In particular two paths of such ISPs are chosen in the Northern 
gyre in the time intervals [0,339] and [446,880]. From each of these paths, a 
DHT is computed which is in the same geographical area although its coor- 
dinates are determined for a different time range. In Fig. [20] we show the x 
and y evolutions for these trajectories. These coordinates have been computed 
with a different algorithm to that proposed in [T?. Instead each corresponds 
to a trajectory which is in the intersection of a piece of a stable manifold and 
a piece of an unstable manifold which are evolved in backwards and forwards 
time respectively. In this procedure, in order to avoid the numerous intersec- 
tions between stable and unstable manifolds, which make difficult the tracking 
of the trajectory which is distinguished, manifolds are trimmed at each time 
step following the ideas in [11 where a method is described to compute a piece 
of single branch of the stable or of the unstable manifold. This method takes 
advantage of the fact that a DHT must be in the intersection of both manifolds 
at all times, as it is a trajectory, however does not improve the method explained 
in [12] in the sense that it does not allow either to extend the computation of 
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the DHT beyond the time interval in which the ISP exists. Many questions have 
been raised for these trajectories as has been discussed in [TTJIT^. For instance 
as they have been computed only in finite time intervals on which the ISP ex- 
ists, one can ask how to pursue its computation beyond that interval. Another 



open issue in [12 concerns deciding if the two DHT in Fig. 20 computed at 
different times are part of the same trajectory. In |llj . the question is raised 
of whether it can happen that a DHT ceases to be distinguished or hyperbolic. 
In this section we apply our algorithm to compute limit coordinates and verify 
whether trajectories in Fig. [20] are distinguished or not following our definition 
[3] Also we will describe how this definition helps address the questions raised 
in |lll 112] . We have applied our algorithm to compute limit coordinates in the 
domain in which the DHT shown in Fig. |20^ ) exists. In particular we have 
applied it with the input: V = (55, 75) x (1325, 1375) = 120, Im = 300, 

At = 5 days, tq = 2 days. At = 5 days, and 5 = 10^^ km. The time step of 



the Runge-Kutta method is 0.1 days. Fig. 21 1) indicates with a solid line the 



projection onto the x ^ y plane of the trajectory depicted in figure 20 1) in the 



interval (120,300), and with circles the path of limit coordinates. Fig. 21 d) 



shows the evolution of the distances between these trajectories. This confirms 



that the trajectory displayed in Fig. 20 1) is also distinguished in the sense of 



definition [3] in the time interval [120, 330] with accuracy e = 8 - 10 ^ km. Thus in 
this time interval, limit coordinates give a method for computing DT different 



from those proposed in [T^ [T7]. Circles in Fig. 22 show the location versus 
time of the x limit coordinates computed with our algorithm. The solid line 
represents a trajectory obtained after integrating with a 5th order Runge-Kutta 
method forwards and backwards in time the initial condition of the circle at day 
285. The dashed line represents the same, but with the initial condition slightly 
perturbed. It is evident that in both cases the trajectories are aligned with 
the path of limit coordinates. The distinguished trajectory is highly hyperbolic 
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backwards in time as in that direction a small perturbation amplifies greatly, 
while it does not do so forwards in time, suggesting that it has a non- hyperbolic 
type of stability in that direction (see comments to Figs. |6]and 15 1. 



Beyond day 300 it is possible to continue the path of limit coordinates. 
Figure |23] shows a diagram at day 330; there is shown the convergence of the 
X component of the minimum of M versus r. This type of convergent diagram 
is not found in this neighbourhood for day 337. On the other hand, although 
it is possible to continue the path of limit coordinates beyond day 300, Fig. [24| 
proves that this path is not a trajectory. There can be seen the existence of 
different trajectories crossing the path, confirming that it is not a trajectory as 
otherwise it would violate the uniqueness of the solution. Therefore, following 
our construction it is possible to say that beyond day 300 the trajectory is no 
longer distinguished. 



Fig. 25 confirms that the trajectory in Fig. [20p) is also distinguished in 
the sense of definition |3] in the time interval (470, 860) with accuracy e = 3 km. 
In particular to compute the path in Fig. [25] we have applied the algorithm of 
section 3 with the input: V ^ (50, 65) x (1255, 1270) km^, to = 470, = 860, 
At = 5 days, tq — 2 days. At = 7 days, and 5 = 10^^ km. The Runge-Kutta 
time step is 0.1 days. In the time interval from day 600 to day 650 some of the 
input parameters were modified as follows: V — (73.5, 75.5) x (1384, 1392) km'^. 
To — 40 days, and At — 1 day. This was due to the presence of nearby elliptic 
type minima in the function M , that made difficult tracking the path of the 
limit coordinates with the previous input. 

Finally, we discuss the existence of non-hyperbolic distinguished trajectories 
in this data set. The presence of this type of trajectories has not been addressed 
before, and we do not have any benchmark solution. We have looked for this 
type of trajectory in areas of the flow where Eulerian eddies seemed to persist 
for long times. Figure [26] represents the function M at day 370 for t = 150 
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and T = 250. In these figures there can be seen the structure of an eddie at 



the centre even for rather long r- values, however figure 27 does not confirm 
the convergence of the minimum of M towards a constant value. On the other 
hand, the slow convergence in diagram[T3]towards the non-hyperbolic trajectory, 
already suggested that long time intervals were required for that purpose, and 
those intervals might be difficult to find in realistic fiows such like the one 
analysed here, in which one is provided just with a finite time datat set. 

6 Conclusions 

In this paper we have proposed a new definition of distinguished trajectory that 
attempts to extend the concept of fixed point and periodic orbit to aperiodic 
dynamical systems. The concept of fixed point is trivially contained in the def- 
inition. Regarding other especially useful trajectories in dynamical systems, for 
instance periodic orbits, we have not proved that they fall within the definition 
in a general way, but we have numerically verified it for selected 2D and 3D 
examples. The definition can be implemented numerically and the procedure 
consists in determining a path of limit coordinates. We have analysed exact 
examples for the Duffing equation with known distinguished trajectories, both 
periodic an aperiodic, and we have found that the path of limit coordinates 
coincides, to within numerical accuracy, with the distinguished trajectories and 
therefore those trajectories are identified also as distinguished in the framework 
of our definition. Our definition is novel with respect to previous works deal- 
ing with distinguished trajectories, because it is applicable to non- hyperbolic 
trajectories. In particular we have studied a periodic orbit of the Duffing equa- 
tion with non- hyperbolic stability and is also recognised as distinguished by our 
definition. In this case the function AI from which the limit coordinates are 
computed seems to be a suggestive tool for characterising Lagrangian eddies. 
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We have tested our definition in the context of reaHstic aperiodic flows where 
distinguished hyperbohc trajectories had been found [TTJ [T^ . Again we have 
identified these trajectories by paths of hmit coordinates in certain time inter- 
vals. Beyond these time intervals the trajectories are no longer distinguished 
according to our definition. Thus in the context of the definitions provided in 
this paper, the property of a trajectory of being distinguished may be lost in 
time. Also we have found evidence that the hyperbolicity of these trajectories 
is not constant in time. These two statements provide answers to the open 
questions mentioned in the text that have been addressed in [TTJ [T2] . 

Acknowledgements 

The authors have been supported by CSIC grants PI-200650I224 and OCEAN- 
TECH (No. PIF06-059), Consohder I- MATH (C3-0104) and the Comunidad de 
Madrid project SIMUMAT S-0505-ESP-0158. It is a pleasure to acknowledge 
many useful conversations with Steve Wiggins, Des Small, Peter Haynes, Emilio 
Hernandez-Garcfa, Cristobal Lopez, Antonio Turiel, Emilio Garcia-Ladona, Michal 
Branicki, Wenbo Tang, J.J. L. Velazquez on numerous issues related to this 
project. We also thank very valuable suggestions from Daniel Fox, Andrew 
Thompson, and Jodie Holdway. The computational part of this work was 
done using the CESGA computers SVGD and FINIS TERRAE and using the 
SIMUMAT-CSIC cluster ODISEA. 

References 

[1] Aref, H., 1984. Stirring by chaotic advection. J. Fluid Mech. 143, 1 21. 

[2] Khakhar DV, Ottino J.M. Fluid mixing (stretching) by time-periodic se- 
quences for weak flows. Physics of Fluids 29 11 3503-3505 (1986). 



29 



[3] Dombre T., Prisch U., Greene J.M., Henon M., Mehr A., Soward A.M. 
Chaotic Streamlines in the ABC flows Journal of Fluid Mechanics 167, 353- 
391 (1986). 

[4] Wiggins, S.: Chaotic Transport in Dynamical Systems, Springer- Verlag, 
New York, 1992 

[5] Wiggins, S.: Introduction to Applied Nonlinear Dynamical Systems and 

Chaos , Springer- Verlag, New York, 2003. 

[6] Guckenheimer, J., Holmes, P.: Nonlinear Oscillations, Dynamical Systems, 
and Bifurcations of Vector Fields, Springer- Verlag, New York, 2002 

[7] Malhotra, N. and Wiggins, S.: Geometric Structures, Lobe Dynamics, and 
Lagrangian Transport in Flows with Aperiodic Time Dependence, with Ap- 
pUcations to Rossby Wave Flow, J. Nonlinear Sci., 8, 401456, 1998. 

[8] Haller, G. and Poje, A.: Finite time transport in aperiodic flows, Physica 
D, 119, 3/4, 352380, 1998. 

[9] Langa JA, Robinson JC, Suarez A. Stability, instability, and bifurcation 
phenomena in non-autonomous differential equations. Nonlinearity 15 (3) 887- 
903 (2002). 

[10] Langa JA, Robinson JC, Suarez A. Bifurcations in non-autonomous scalar 
equations. Journal of Differential Equations. 221 (1) 1-35 (2006). 

[11] Mancho A. M., Small D., Wiggins S., A tutorial on dynamical systems 
concepts applied to Lagrangian transport in oceanic flows dcflncd as flnite 
time data sets: Theoretical and computational issues. Physics Reports 437, 
3-4 (2006). 



30 



[12] Mancho, A. M., Small, D., Wiggins, S., 2004. Computation of hyperbolic 
and their stable and unstable manifolds for oceanographic flows represented 
as data sets. Nonlinear Process. Geophys. 11, 17 33. 

[13] Mancho A. M., Hernandez-Garcia E., Small D., Wiggins S., Fernandez V. 
Lagrangian transport through an ocean front in the North- Western Mediter- 
ranean Sea. Journal of Physical Oceanography 38, 6, 1222-1237 (2008). 

[14] AurcU E., Boffctta G., Crisauti A., Paladin G., Vulpiani A. Predictability 
in the large: an extension of the concept of Lyapunov exponent. J. Phys. A; 
Math. General, 30: 1-26, 1997. 

[15] Haller, G. Distinguished Material surfaces and coherent structures in three- 
dimensional fluid flows. Physica D, 149: 248-277, 2001. 

[16] Nese JM, Quantifying local predictability in phase space. Physica D 35: 
237-250, 1989. 

[17] Ide, K., Small. D., Wiggins, S., 2002. Distinguished hyperbolic trajectories 
in time dependent fluid flows: analytical and computational approach for 
velocity flelds deflned as data sets. Nonlinear Process. Geophys. 9, 237 263. 

[18] Ju, N., Small, D., Wiggins, S., 2003. Existence and computation of hy- 
perbolic trajectories of aperiodically time-dependent vector flelds and their 
approximations. Int. J. Bif. Chaos 13, 1449 1457. 

[19] G. Haller. Exact theory of unsteady separation for two dimensional flows. 
Journal of Fluid Mechanics 512, 257-311 (2004) 

[20] Eisenbach S., Friedrich R., Large-eddy simulation of flow separation on an 
airfoil at a high angle of attack and Re=10(5) using Cartesian grids. Theo- 
retical and Computational Fluid Dynamics 22 (3-4) 213-225 (2008). 



31 



[21] Coulliette, C, Wiggins, S., 2001. Intergyre transport in a wind-driven, 
quasigeostrophic double gyre: an application of lobe dynamics. Nonlinear 
Process. Geophys. 8, 69 94. 

[22] Szeri, A., Leal, L. G., and Wiggins, S.: On the Dynamics of Suspended Mi- 
crostructure in Unsteady, Spatially Inhomogeneous Two-Dimensional Fluid 
Flows, J. Fluid Mech., 228, 207241, 1991. 

[23] Mancho, A. M., Small, D., Wiggins, S., Idc, K.,2003. Computation of Sta- 
ble and Unstable Manifolds of Hyperbolic Trajectories in Two-Dimensional, 
Aperiodically Time-Dependent Vectors Fields. Physica D 182, 188-222. 

[24] Mancho, A. M., Small, D., Wiggins, S., Ide, K.,2006. A comparison of meth- 
ods for interpolating chaotic flows from discrete velocity data. Computers and 
Fluids 35, 416-428. 

[25] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical 
Recipes in Fortran 77, The Art of Scientiflc Computing, 2nd ed., Cambridge 
University Press, Cambridge, 1999. 

[26] Dritschel, D.G., Contour dynamics and contour surgery: numerical algo- 
rithms for extended, high-resolution modelling of vortex dynamics in two- 
dimensional, inviscid, incompressible flows, Comput. Phys. Rep. 10 (1989) 
77146. 

[27] Branicki, M, Mancho, A. M., Wiggins, S. A Lagrangian description of trans- 
port associated with a Front-Eddy interaction: application to data from the 
North- Western Mediterranean Sea. Preprint submitted for publication. 

[28] D. B. Chelton, M. G. Schlax, R. M. Samelson, and R. A. de Szoeke, Global 
observations of large oceanic eddies. Geophysical Research Letters 34 (2007), 
L15606. 



32 



[29] Haller G, An objective definition of a vortex, Journal of Fluid Mechanics 
525, 1-26 (2005). 



33 



Figure captions 

Figure 1 

Solutions ^ for different initial conditions x* . a) Solutions for positive 
times t > 0. b) Solutions for positive and negative times. 
Figure 2 

a) Function M{x*)t=o_r evaluated over the solutions ([6]). Dashed line t = 3, 
solid line t = 4. b) Position of the x*-coordinate at the minimum of the function 
Mt=o,T as a function of r. The horizontal dashed line marks the position of the 
DHT. 

Figure 3 

Position of the -coordinate at the minimum of the function Aft=o,r as a 
fimction of t. The function Mf=o,r is considered for the solutions in (|8|. The 
horizontal dashed line marks the position of the DHT. 

Figure 4 

Contour plot of the function A/(^o.t=2(x) in the open set V = (—0.2, 0.2) x 
(—0.2, 0.2). The minimum corresponds to the black tone. 
Figure 5 

Contour plot of the function M in the open set x S (—0.2, 0.2) x (—0.2, 0.2). 
a) Mi^o,r=5(x); b) Mt=o,r=io(x); c) M4=o,r=5o(x). 
Figure 6 

a) x-coordinate versus time for the DHT (thick solid line) and those trajec- 
tories integrated with initial conditions at the global minima of Figs. [5^) (solid 
line) and b) (dashed line); b) x-coordinate versus time for the DHT (thick solid 
line), a trajectory integrated with initial condition at the global minimum of 
Fig. ^) (solid line) and a trajectory integrated at a non-global but relative 
minimum of the same figure (dashed line). 

Figure 7 
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Evolution of the coordinates of the global minimum of M versus t. a) The x 
coordinate; b) the y coordinate. These plots show the convergence to the DHT 
whose position is marked with a dashed horizontal line. 

Figure 8 

A (5-grid in M^. The centre or interior point is marked with the white dot. 
Figure 9 

A schematic projection onto the plane of a possible trajectory from — r 
to r with L points. 
Figure 10 



a) Representation of both the distinguished hyperbolic trajectory (11) and 
its approximation obtained with the proposed numerical algorithm for e = 0.1; 
b) distance between the exact and the numerical approach. 

Figure 11 

Representation of the maximum r required to approach the DT to within 
accuracy 5 — 10^^ versus time. 
Figure 12 

Contour plot of the function M in the open set x e (—1.2, — 0.8) x (—0.2, 0.2)]. 
a) M4^o,T=io(x); b) Af4^o,T=3oo(x). 
Figure 13 

Evolution of the coordinates of the global minimum of AI versus t aX — 0. 
a) The x coordinate; b) the y coordinate. These plots show the convergence 
of the minima to the coordinates of the DET whose position is marked with a 
continuous horizontal line. 

Figure 14 

a) Dotted line represents the exact non-hyperbolic distinguished trajectory 
and the solid line stands for the numerically computed limit coordinates; b) 



distance between the exact non- hyperbolic trajectory (20 1 and the limit coordi- 
nates. 
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Figure 15 

a) x-coordinate versus time for the DET (solid line) and the trajectory inte- 
grated taking as initial data the limit coordinates located at time to — (dashed 
line); b) time evolution of the differences between these trajectories. 

Figure 16 

a) Dashed line represents the exact distinguished hyperbolic trajectory of the 
rotating Duffing equation and the solid line stands for the numerically computed 
one; b) distance between the exact and the numerical distinguished hyperbolic 
trajectories. 

Figure 17 

Evolution of the coordinates of the global minimum of M for the 3D example 
versus r at io = 0- The x coordinate; b) the y coordinate; c) the z coordinate. 
These plots show the convergence of the minima to the coordinates of the DHT 
whose position is marked with a dashed horizontal line. 

Figure 18 

a) The solid line represents the exact distinguished hyperbolic trajectory 
of the 3D equation and circles stand for numerically computed coordinates; b) 
distance between the exact and the numerical distinguished hyperbolic trajec- 
tories. 

Figure 19 

Contour plot of the streamfunction produced by the quasigeostrophic model 
at day 300. 
Figure 20 

Distinguished hyperbolic trajectories in the Northern gyre of the quasi- 
geostrophic model reported in [12 . a) Evolution of the x and y coordinates 
in the time interval [5, 338]; b) evolution of the x and y coordinates in the time 
interval [450, 880] 

Figure 21 
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a) Solid line represents the projection on the phase space of the distinguished 



hyperbolic trajectory depicted in Fig. 20 1) and circles stands for the numerically 



computed limit coordinates; b) distance between the trajectories represented in 
a). 

Figure 22 

Circles stand for the x component of the limit coordinates in the time range 
where they approach a DT. The solid line represents a trajectory integrated 
with a 5th order Runge-Kutta method passing the through limit coordinates at 
day 285. The dashed line is a trajectory integrated from the same condition 
plus a small perturbation. 

Figure 23 

a) X component of the minimum of M versus r at day 330; b) y component 
of the minimum of M versus r at the same day. 
Figure 24 

a) Circles stand for the x component of the limit coordinates versus time 
and the solid lines stand for different trajectories; b) the same as a) but for the 
y component. 

Figure 25 

a) Solid line represents the projection on the phase space of the distinguished 



hyperbolic trajectory depicted in Fig. 20 and circles stands for the numerically 



computed limit coordinates; b) distance between the trajectories represented in 
a). 

Figure 26 

a) Contour plot of Mi=37o,r=i50i the elliptic minimum is in the dark area 
almost at the centre; b) contour plot of Aft=37o,r=250- 
Figure 27 

a) X component of the minimum of M versus r at day 370; b) y component 
of the minimum of M versus r at the same day. 
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Table captions 

Table 1 

Relative errors for several ellipse lengths, computed with a linear interpola- 
tion over L points on the curve. 
Table 2 

Relative errors for several ellipse lengths, computed with Dritschel interpo- 
lation over L points on the curve. 
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L 


Linear interpolation 


Ratio between axes 


1 


2 


5 


10 


100 


1000 


10 


8.16 


8.80 


9.48 


9.73 


9.99 


10.00 


10^ 


2.63 


3.36 


3.86 


3.99 


4.05 


4.05 


10-^ 


0.83 


1.08 


1.24 


1.29 


1.31 


1.31 


10* 


0.26 


0.34 


0.39 


0.41 


0.41 


0.41 


10*= 


8.34 X 10-^ 


0.11 


0.12 


0.13 


0.13 


0.13 


10** 


2.64 X 10"^ 


3.42 X IQ-^ 


3.94 X 10"^ 


4.08 X 10"^ 


4.14 X 10"^ 


4.14 X 10"^ 


10' 


8.34 X IQ-^ 


1.08 X 10"^ 


1.25 X 10-^ 


1.29 X 10"^ 


1.31 X 10"^ 


1.31 X 10"^ 


10« 


2.64 X IQ-'' 


3.42 X IQ-^ 


3.95 X 10"^ 


4.08 X 10"^ 


4.14 X 10"^ 


4.14 X 10"^ 


10« 


8.34 X 10-" 


1.08 X 10-^ 


1.25 X 10-^ 


1.29 X 10-^ 


1.31 X 10-^ 


1.31 X 10-^ 
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L 


Dritschel interpolation 


Ratio between axes 


1 


2 


5 


10 


100 


1000 


10 


0.99 


0.67 


0.33 


0.25 


0.26 


0.27 


10^ 


3.49 X 10-2 


1.37 X 10-2 


4.47 X 10--^ 


2.81 X 10-^ 


2.08 X 10-'^ 


2.61 X 10--* 


10^ 


1.11 X 10"''' 


3.83 X 10"* 


9.20 X 10-^ 


4.43 X 10-" 


2.23 X 10--' 


1.98 X 10-" 


10* 


3.53 X 10-*= 


1.17 X 10"" 


4.44 X 10"" 


5.23 X 10"" 


2.80 X 10" ' 


7.01 X 10"' 


10" 


1.12 X 10-" 


3.64 X 10" ' 


2.17 X 10"" 


4.48 X 10-" 


5.45 X 10"** 


9.20 X 10" ' 


10" 


3.51 X 10-*^ 


1.14 X 10-*^ 


2.10 X 10-" 


4.46 X 10-" 


5.222 X 10-** 


9.22 X 10-'' 


10'^ 


1.11 X lO"'* 


3.68 X 10-"' 


2.10 X 10-'' 


4.46 X 10-" 


5.21 X 10-** 


9.22 X 10- ' 


10« 


2.18 X 10"^^ 


1.30 X 10"^^ 


2.10 X 10-" 


4.46 X 10-" 


5.22 X 10"** 


9.22 X 10" ' 


10^ 


3.72 X 10"^" 


3.37 X 10"^" 


2.10 X 10"" 


4.46 X 10"" 


5.22 X 10"** 


9.25 X 10" ' 



Table 2: 
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